Read and wrangle in data

output_main <- read.csv(here("data","bivalve_model","LookUpTable_201416.csv")) %>% 
  dplyr::select("site", "long", "lat", "yield")
output_s1 <- read.csv(here("data","bivalve_model","LookUpTable_s1.csv"))
output_s2 <- read.csv(here("data","bivalve_model","LookUpTable_s2.csv"))
output_s3 <- read.csv(here("data","bivalve_model","LookUpTable_s3.csv"))
s1 <- output_s1 %>% 
  mutate(s1_yield = yield) %>% 
  dplyr::select("site", "long", "lat", "s1_yield")

s2 <- output_s2 %>% 
  mutate(s2_yield = yield) %>% 
  dplyr::select("site", "long", "lat", "s2_yield")

s3 <- output_s3 %>% 
  mutate(s3_yield = yield) %>% 
  dplyr::select("site", "long", "lat", "s3_yield")

Merging the 3 frames. Format will be wide so then we will make it longer.

wide1 <- merge(output_main, s1, by = c("site", "long", "lat"), all = TRUE)

wide2 <- merge(s2, s3, by = c("site", "long", "lat"), all = TRUE)

wide <- merge(wide1, wide2, by = c("site", "long", "lat"), all = TRUE)
alldata <- wide %>% 
  pivot_longer(yield:s3_yield,
               names_to = "test",
               values_to = "yield")

Visualizations

##boxplot and violin plot
ggplot(data = alldata, aes(x = test, y = yield, fill = test)) +
  geom_boxplot() +
  #geom_jitter(color="black", size=0.4, alpha=0.9) +
  theme(
      legend.position="none",
      plot.title = element_text(size=11)
    ) +
  theme_minimal()

ggplot(data = alldata, aes(x = test, y = yield, fill = test)) +
  geom_violin() +
  #geom_jitter(color="black", size=0.4, alpha=0.9) +
  # theme(
  #     legend.position="none",
  #     plot.title = element_text(size=11)
  #   ) +
  # theme(legend.title = element_blank()) +
  xlab("") +
  ylab("Yield (lbs)") +
  theme_minimal() +
  theme(legend.position = "none")

#ggsave(filename = here("figures", "aquaculture_sensitivity.jpeg"), width = 10, height = 8, units = "in")

Percent change

wide_percents <- wide %>% 
  mutate(s1_percchange = (((s1_yield-yield)/yield)*100),
         s2_percchange = (((s2_yield-yield)/yield)*100),
         s3_percchange = (((s3_yield-yield)/yield)*100),
         )

summary_percents <- wide_percents %>% 
  #group_by(site) %>% 
  summarize(mean_s1percchhange = mean(s1_percchange),
            mean_s2percchhange = mean(s2_percchange),
            mean_s3percchhange = mean(s3_percchange)
            )

ANOVA

summary(alldata)
##       site           long             lat            test          
##  Min.   :   6   Min.   :-124.6   Min.   :32.35   Length:18396      
##  1st Qu.:1240   1st Qu.:-122.8   1st Qu.:33.26   Class :character  
##  Median :2415   Median :-119.7   Median :33.93   Mode  :character  
##  Mean   :2476   Mean   :-120.5   Mean   :35.33                     
##  3rd Qu.:3757   3rd Qu.:-118.7   3rd Qu.:37.67                     
##  Max.   :5000   Max.   :-117.1   Max.   :42.01                     
##      yield         
##  Min.   :  871687  
##  1st Qu.: 2836668  
##  Median : 3634054  
##  Mean   : 4433778  
##  3rd Qu.: 5955865  
##  Max.   :11570374
aov <- aov(yield ~ test, 
           data = alldata)

summary(aov)
##                Df    Sum Sq   Mean Sq F value Pr(>F)    
## test            3 2.520e+16 8.400e+15    2831 <2e-16 ***
## Residuals   18392 5.458e+16 2.968e+12                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hist(aov$residuals)

ggplot(data = alldata, aes(sample = yield)) +
  geom_qq() +
  facet_wrap( ~ test) +
  theme_minimal()

Visualization data transformations

This section deals with the production data and transformed production data from each of the 4 sub-models. We want to compare the transformed data to the theoretical transformation (line resulting from the equation we used to transform the data - how well do the data points sit on this line?).

I also want to use this opportunity to look at the production data (pre-transformation) in some non-spatial ways.

I will read in the rasters and then extract the value in each cell.

Read in data using terra

library(terra)
library(sf)

Wave energy

wave <- terra::rast(here("data", "tiff", "wave_linfunc_v1.tif"))
wave_raw <- terra::rast(here("data", "tiff", "wave_notransform.tif"))

terra::plot(wave, main = "Wave Energy (transformed linearly)")

terra::plot(wave_raw, main = "Wave Energy (no transformation)")

View head and tail of data frame

wave_df <- terra::as.data.frame(wave, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)

head(wave_df)
##              x        y    Band_1
## 4884 -124.5397 41.97784 0.8077396
## 4886 -124.4588 41.97784 0.7711606
## 4887 -124.4184 41.97784 0.7387354
## 4888 -124.3780 41.97784 0.6629785
## 4889 -124.3376 41.97784 0.6135752
## 4890 -124.2971 41.97784 0.5418847
tail(wave_df)
##               x        y    Band_1
## 51040 -119.0416 32.43699 0.6873534
## 51041 -119.0012 32.43699 0.6822757
## 51235 -119.0416 32.39656 0.6925157
## 51236 -119.0012 32.39656 0.6886186
## 51237 -118.9607 32.39656 0.6806911
## 52815 -118.2330 32.07314 0.5573189
ggplot(data = wave_df, aes(x = Band_1)) +
  geom_histogram() +
  theme_minimal()

###################################################

wave_raw_df <- terra::as.data.frame(wave_raw, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)

head(wave_raw_df)
##              x        y wave_notransform
## 4884 -124.5397 41.96992         38814.79
## 4886 -124.4588 41.96992         36715.67
## 4887 -124.4184 41.96992         34854.92
## 4888 -124.3780 41.96992         30507.54
## 4889 -124.3376 41.96992         27672.48
## 4890 -124.2971 41.96992         23558.46
tail(wave_raw_df)
##               x        y wave_notransform
## 51040 -119.0416 32.43665         31906.31
## 51041 -119.0012 32.43665         31614.92
## 51235 -119.0416 32.39625         32202.55
## 51236 -119.0012 32.39625         31978.92
## 51237 -118.9607 32.39625         31523.99
## 52815 -118.2330 32.07309         24444.17
ggplot(data = wave_raw_df, aes(x = wave_notransform)) +
  geom_histogram() +
  theme_minimal()

# Creating a scatterplot of data with the transformation line on it

graph_wa <- ggplot(data = wave_df, aes(x = , y =)) +
  geom_point() +
  theme_minimal()

Merging the pre-transformed and transformed data frames

####### YIKES THIS DOESNT WORK ################
# wave_all <- merge(wave_raw_df, wave_df, by = c("x", "y"), all = FALSE)

#trying to combine spatraters first and then will try to convert to tabular data from this format

wave_sprc <- sprc(wave_raw, wave)
# terra::plot(wave_sprc, main = "Testing sprc")

wave_sprc_df <- terra::as.data.frame(wave_sprc, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)

Wind energy

wind <- terra::rast(here("data", "tiff", "wind_linfunc_v_ProjectRas.tif"))
wind_raw <- terra::rast(here("data", "tiff", "wind_notransform.tif"))

terra::plot(wind, main = "Wind Energy (transformed linearly)")

terra::plot(wind_raw, main = "Wind Energy (no transformation)")

View head and tail of data frame

wind_df <- terra::as.data.frame(wind, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)

head(wind_df)
##                 x        y    Band_1
## 1083981 -124.3807 42.06418 0.7324007
## 1083984 -124.3708 42.06418 0.7286220
## 1083985 -124.3675 42.06418 0.7284269
## 1083988 -124.3577 42.06418 0.7237991
## 1083989 -124.3544 42.06418 0.7185299
## 1083991 -124.3478 42.06418 0.7152024
tail(wind_df)
##                  x        y    Band_1
## 10049327 -119.0085 32.16150 0.4130057
## 10052294 -119.0381 32.15821 0.4195321
## 10052295 -119.0348 32.15821 0.4183595
## 10052296 -119.0315 32.15821 0.4169410
## 10052297 -119.0282 32.15821 0.4160608
## 10052298 -119.0250 32.15821 0.4160352
ggplot(data = wind_df, aes(x = Band_1)) +
  geom_histogram() +
  theme_minimal()

###################################################

wind_raw_df <- terra::as.data.frame(wind_raw, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)

head(wind_raw_df)
##                x        y   Band_1
## 990023 -124.3441 42.06536 23472.07
## 990026 -124.3342 42.06536 23263.02
## 990027 -124.3309 42.06536 23065.66
## 990029 -124.3243 42.06536 23032.67
## 990030 -124.3210 42.06536 22990.36
## 990031 -124.3177 42.06536 22983.74
tail(wind_raw_df)
##                 x        y   Band_1
## 9913185 -119.0245 32.16039 15022.49
## 9913186 -119.0212 32.16039 15014.09
## 9913187 -119.0179 32.16039 15005.87
## 9913188 -119.0146 32.16039 14989.09
## 9916144 -119.0343 32.15710 15081.72
## 9916145 -119.0310 32.15710 15048.69
ggplot(data = wind_raw_df, aes(x = Band_1)) +
  geom_histogram() +
  theme_minimal()

Merging

# wind_all <- merge(wind_raw_df, wind_df, by = c("x", "y"), all = FALSE)

Aquaculture

aqua <- terra::rast(here("data", "tiff", "aqua_sfunc_v2.tif")) #???????
aqua_raw <- terra::rast(here("data", "tiff", "aqua_notransform.tif"))

terra::plot(aqua, main = "Aquaculture Production (S-transformation)")

terra::plot(aqua_raw, main = "Aquaculture Production (no transformation)")

View head and tail of data frame

aqua_df <- terra::as.data.frame(aqua, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)

head(aqua_df)
##              x        y    Band_1
## 8693 -124.4787 42.00743 0.5604585
## 8694 -124.4487 42.00743 0.6017970
## 8695 -124.4187 42.00743 0.6330283
## 8696 -124.3888 42.00743 0.6407387
## 8697 -124.3588 42.00743 0.6946170
## 8698 -124.3288 42.00743 0.7214648
tail(aqua_df)
##               x        y       Band_1
## 93034 -119.0533 32.41554 4.251122e-04
## 93035 -119.0233 32.41554 8.228081e-05
## 93036 -118.9933 32.41554 9.777889e-06
## 93297 -119.0533 32.38557 1.366302e-04
## 93299 -118.9933 32.38557 2.556847e-06
## 93300 -118.9633 32.38557 1.362262e-06
ggplot(data = aqua_df, aes(x = Band_1)) +
  geom_histogram() +
  theme_minimal()

###################################################
aqua_raw_df <- terra::as.data.frame(aqua_raw, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)

head(aqua_raw_df)
##              x        y aqua_notransform
## 8693 -124.4787 41.99691          7521588
## 8694 -124.4487 41.99691          7716453
## 8695 -124.4187 41.99691          7870482
## 8696 -124.3888 41.99691          7909508
## 8697 -124.3588 41.99691          8194795
## 8698 -124.3288 41.99691          8346378
tail(aqua_raw_df)
##               x        y aqua_notransform
## 93034 -119.0533 32.41518          3064938
## 93035 -119.0233 32.41518          2994500
## 93036 -118.9933 32.41518          2958244
## 93297 -119.0533 32.38523          3010470
## 93299 -118.9933 32.38523          2948923
## 93300 -118.9633 32.38523          2946289
ggplot(data = aqua_raw_df, aes(x = aqua_notransform)) +
  geom_histogram() +
  theme_minimal()

Merging

aqua_all <- merge(aqua_raw_df, aqua_df, by = c("x", "y"), all = FALSE)

Fisheries

SOMETHNG IS STILL WRONG HERE. I think maybe there might be an issue with the pre-transformation raster still. Transformed data has 1,640,319 observations and pre-transformation only has 15,541, which is weird.

fishy <- terra::rast(here("data", "tiff", "catchdata_outp_ProjectRas.tif"))
fishy_raw <- terra::rast(here("data", "tiff", "catchdata_notransform.tif"))

terra::plot(fishy, main = "Fisheries Production (Z-transformation)")

terra::plot(fishy_raw, main = "Fisheries Production (no transformation)")

View head and tail of data frame

fishy_df <- terra::as.data.frame(fishy, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)

head(fishy_df)
##                x        y catchdata_outp_ProjectRas
## 940787 -125.5002 41.99897                 0.6294906
## 940788 -125.4969 41.99897                 0.6294906
## 940789 -125.4936 41.99897                 0.6294906
## 940790 -125.4903 41.99897                 0.6294906
## 940791 -125.4869 41.99897                 0.6294906
## 940792 -125.4836 41.99897                 0.6294906
tail(fishy_df)
##                 x        y catchdata_outp_ProjectRas
## 9273852 -117.1849 32.33468                 0.8797805
## 9273853 -117.1816 32.33468                 0.8797805
## 9273854 -117.1782 32.33468                 0.8797805
## 9273855 -117.1749 32.33468                 0.8797805
## 9273856 -117.1716 32.33468                 0.8797805
## 9273857 -117.1683 32.33468                 0.8797805
ggplot(data = fishy_df, aes(x = catchdata_outp_ProjectRas)) +
  geom_histogram() +
  theme_minimal()

###################################################
fishy_raw_df <- terra::as.data.frame(fishy_raw, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)

head(fishy_raw_df)
##           x        y catchdata_notransform
## 7  -125.479 41.99396              25653.79
## 8  -125.445 41.99396              25653.79
## 9  -125.411 41.99396              25653.79
## 10 -125.377 41.99396              25653.79
## 11 -125.343 41.99396              25653.79
## 12 -125.309 41.99396              25653.79
tail(fishy_raw_df)
##              x        y catchdata_notransform
## 72054 -118.849 32.33796            211535.130
## 72099 -117.319 32.33796              2010.012
## 72100 -117.285 32.33796              2010.012
## 72101 -117.251 32.33796              2010.012
## 72102 -117.217 32.33796              2010.012
## 72103 -117.183 32.33796              2010.012
ggplot(data = fishy_raw_df, aes(x = catchdata_notransform)) +
  geom_histogram() +
  theme_minimal()

 # xlim(c(-1, 500000))
fishy_all <- merge(fishy_raw_df, fishy_df, by = c("x", "y"), all = FALSE)

Co-location Suitability Model

colo <- terra::rast(here("data", "tiff", "CellSta_colocation_off.tif"))

terra::plot(colo, main = "Co-location Suitability Score")

colo_df <- terra::as.data.frame(colo, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)

head(colo_df)
##             x        y    Band_1
## 490 -125.4974 41.98627 0.6294906
## 491 -125.4570 41.98627 0.7640004
## 492 -125.4165 41.98627 0.7640711
## 493 -125.3761 41.98627 0.7652843
## 494 -125.3357 41.98627 0.7650780
## 495 -125.2952 41.98627 0.7662386
tail(colo_df)
##               x        y    Band_1
## 56103 -119.0694 32.36456 0.3808336
## 56104 -119.0290 32.36456 0.3760225
## 56105 -118.9886 32.36456 0.3719020
## 56106 -118.9482 32.36456 0.3670245
## 56107 -118.9077 32.36456 0.3618057
## 56108 -118.8673 32.36456 0.3588376
colo_hist <- ggplot(data = colo_df, aes(x = Band_1)) +
  geom_histogram() +
  theme_minimal() +
  xlab("Co-location Suitability Score") +
  ylab("Count")
# ggsave(filename = here("figures", "colocation_histogram.jpeg"), width = 8, height = 8, units = "in")

Scatter plot showing the co-location suitability score related to latitude

# library(hrbrthemes) #having issues downloading some of these dependencies

ggplot(data = colo_df, aes(x = y, y = Band_1)) +
  geom_point(size = 1) +
  # geom_smooth(method="lm") +
  theme_minimal() +
  xlab("Latitude") +
  ylab("Co-location Suitability Score")

# ggsave(filename = here("figures", "colocation_latitude.jpeg"), width = 8, height = 8, units = "in")
############################################################################

ggplot(data = colo_df, aes(x = x, y = Band_1)) +
  geom_point(size = 1) +
  # geom_smooth(method="lm") +
  # geom_line(data = long_lm2, aes( y = y)) +
  theme_minimal() +
  xlab("Longitude") +
  ylab("Co-location Suitability Score")

# ggsave(filename = here("figures", "colocation_longitude.jpeg"), width = 8, height = 8, units = "in")

Linear models to correspond each graph above

long_lm <- lm(Band_1 ~ y, data = colo_df)
#long_lm2 <- glm(Band_1 ~ y, family = Gamma(link = "log"), data = colo_df)

lat_lm <- lm(Band_1 ~ x, data = colo_df)

library(broom)
# broom::tidy(long_lm)
# broom::tidy(lat_lm)
glance(long_lm)
## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic p.value    df logLik    AIC    BIC
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl>  <dbl>  <dbl>
## 1     0.192         0.192 0.177     2591.       0     1  3437. -6869. -6847.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>
glance(lat_lm)
## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic p.value    df logLik    AIC    BIC
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl>  <dbl>  <dbl>
## 1     0.385         0.385 0.154     6835.       0     1  4929. -9853. -9831.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>
summary(long_lm)
## 
## Call:
## lm(formula = Band_1 ~ y, data = colo_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50800 -0.10964  0.02131  0.12494  0.49704 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.4783852  0.0215951  -22.15   <2e-16 ***
## y            0.0300670  0.0005907   50.90   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1767 on 10924 degrees of freedom
## Multiple R-squared:  0.1917, Adjusted R-squared:  0.1916 
## F-statistic:  2591 on 1 and 10924 DF,  p-value: < 2.2e-16
summary(lat_lm)
## 
## Call:
## lm(formula = Band_1 ~ x, data = colo_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.51035 -0.10076  0.00483  0.10449  0.51011 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.9849297  0.0798722  -74.93   <2e-16 ***
## x           -0.0540747  0.0006541  -82.68   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1541 on 10924 degrees of freedom
## Multiple R-squared:  0.3849, Adjusted R-squared:  0.3848 
## F-statistic:  6835 on 1 and 10924 DF,  p-value: < 2.2e-16

Latitude seems to have the better correlation and higher R value

plot(lat_lm)

Break point

library(nlme)
library(ciTools)

model_long <- predict.glm(best.mod)

Visualizations

# wave_graph <- ggplot(data = wave_all, aes(x = Band_1, y = wave_notransform)) +
#   geom_point() +
#   theme_minimal()